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ABSTRACT 

A  method  is  described  for  calculating  the  flow  of  an 
invlscld,  Ideal  gas  between  a  blunt  body  moving  at  high 
Mach  number  and  the  detached  shock  that  precedes  It.   The 
method  Is  based  on  power-series  expansion  In  two  space 
variables  about  a  point  on  the  shock.   Special  Unlvac 
routines  for  formal  calculation  with  power  series  have 
been  developed?  they  Include  an  Interpretive  routine  of 
such  a  nature  that  the  partial  differential  equations  that 
determine  the  coefficients  of  the  power  series  can  be 
transcribed  directly  Into  a  pseudo-code.   An  Important 
feature  of  the  method  Is  the  use  of  special  floating-decimal 
subroutines  In  which  the  possible  loss  of  significant  digits 
by  cancellation  is  continually  monitored.   Test  calculations 
that  have  been  made  with  the  method  are  reported.   The 
results  appear  to  indicate  that  the  method  has  considerable 
promise  for  fluld-dynamlcal  problems  of  this  kind,  but 
further  development  In  certain  directions  Is  Indicated. 
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NYO-7973 
DETACHED- SHOCK  CALCULATIONS  BY  POWER  SERIES,  I 
Introduction. 


The  feasibility  of  machine  calculations  by  power 
series  in  a  certain  class  of  fluid-dynamical  problems 
has  been  tested  on  the  Univac.   The  aim  was  to  calculate 
the  steady,  two-dimensional  flow  of  a  compressible  fluid 
when  analytic  Cauchy  data  are  given  on  an  analytic  curve, 
not  a  characteristic.   The  flow-quantities  (pressure,  density, 
velocity  components,  and  the  like)  are  represented  by  power 
series  expansions  in  two  variables,  or  curvilinear  coordinates 
one  of  which  is  constant  on  the  given  curve o   The  point  about 
which  the  expansion  is  made  is  a  point  on  the  curve.   The 
determination  of  the  expansion  coefficients  from  the  given 
data  and  the  partial  differential  equations  of  fluid  dynamics 
is  then  a  formal  calculation,  which  is  made  feasible  by  a 
set  of  subroutines  for  manipulation  of  power  series.   When 
sufficiently  many  coefficients  (say  a  hundred  for  each  flow 
quantity)  have  been  computed,  the  series  are  siimmed  to 
calculate  the  flow  in  a  neighborhood  of  the  point  about 
which  the  expansions  were  made.   The  validity  of  the  ex- 
pansions can  be  tested  empirically  by  l)  comparing  the 
results  obtained  at  a  given  point  from  expansions  about 
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two  or  more  neighboring  points  on  the  curve,  2)  testing 
Bernoulli's  law,  and  J>)    checking  the  constancy  of  entropy 
along  a  streamline. 

The  particular  problem  of  this  class  that  was  chosen 
for  study  Is  to  calculate  the  Invlscld  flow  between  a 
detached  shock  and  a  rotatlonally  symmetric  blunt  body  like 
the  nose  of  a  ballistic  missile  travelling  at  high  Mach 
number  through  air.   The  general  characteristics  of  this 
flow  are  well  known;  the  relationships  between  the  shock, 
the  body,  and  the  streeunllnes,  are  Indicated  In  Figure  1  -- 
a  frame  of  reference  has  been  chosen  In  which  the  body  Is 
at  rest. 


'STREAMLINES 

SHOCK 
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Figure  1 
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The  stand-off  distance  is  of  the  order  of  one  tenth  the 
radius  of  curvature  of  the  body  at  the  nose,  the  actual 
distance  depending  on  the  Mach  number  and  the  assumed 
equation  of  state.   Near  the  axis,  the  flow  is  subsonic; 
farther  away  it  is  supersonic.   Under  usual  conditions, 
the  boundary  layer  that  builds  up  adjacent  to  the  body 
is  of  negligible  thickness;  therefore,  the  assiimption 
of  an  invlscld  flow  throughout  most  of  the  Interesting 
region  is  Justified. 

As  in  several  other  methods  of  attack  on  the  detached- 
shock  problem,  we  start  by  assiiming  a  shape  for  the  shock, 
with  the  Idea  of  determining,  from  the  calculation,  what 
shape  of  body  will  produce  such  a  shock.   If  the  pressure, 
density,  and  velocity,  of  the  air  ahead  of  the  shock  are 
known,  their  values  are  determined  immediately  behind  the 
shock  by  the  Rankine-Hugoniot  Jump  conditions.   If  the 
assumed  shock  surface  and  the  equation  of  state  are  analytic, 
the  flow  quantities  Immediately  behind  the  shock  are  analytic 
and  can  be  expanded  in  powers  of  a  suitable  coordinate  along 
the  shock.   A  preliminary  routine  calculates  the  coefficients 
of  these  expansions;  the  main  routine  calculates  the  re- 
maining coefficients  of  the  full  expansions  in  the  two 
variables  from  the  partial  differential  equations;  a  final 
routine  then  calculates  various  quantities  with  the  aid 
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of  the  power  series  and  prints  them,  for  any  desired  pairs 
of  values  of  the  coordinates. 

An  essential  feature  of  the  calculation  is  the  use 
of  a  special  floating-point  arithmetic,  in  which  there  is 
attached  to  each  nvunber,  in  addition  to  the  fractional  part 
and  exponent,  a  significance  index  which  tells  approximately 
how  many  digits  of  the  fractional  part  are  significant. 
It  appears  that,  owing  to  the  use  of  a  finite  number  of 
decimal  places,  the  series  obtained  by  a  numerical  calcu- 
lation have  some  of  the  characteristics  of  asymptotic 
series:  the  accuracy  improves  for  a  while,  as  the  number  of 
terms  increases,  but  subsequently  deteriorates  rapidly. 
The  explanation  is  that  the  calculation  of  the  coefficients 
is  itself  sufficiently  long  and  complicated  that  loss  of 
significance  can  occur,  with  increasing  severity  for  the 
high-order  coefficients.   The  significance  index  actually 
goes  to  zero  for  sufficiently  high-order  terms.   These 
terms  and,  in  fact,  all  terms  having  indices  below  a 
certain  limit,  must  be  systematically  dropped  from  the 
series  to  obtain  accurate  results.   For  expansions  in  two 
or  more  variables,  in  complicated  problems,  it  would  be 
difficult  to  decide  which  terms  to  retain  without  some 
form  of  significance-monitoring. 

In  Section  2  of  this  report,  the  detached-shock 
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problem  Is  formulated;  In  Section  3,  the  method  of  calcu- 
lation is  described  in  general  terms;  in  section  4,  the 
results  ot   trial  calculations  are  presented.   The  trial 
calculations  were  for  air  (assumed  to  be  an  ideal  gas) 
at  Mach  12,  the  shock  surface  being  taken  as  a  hyperboloid 
of  revolution.   The  floating-point  subroutines  are 
described  in  Appendix  A,  and  the  subroutines  and  pseudo- 
code for  manipulating  power  series  in  Appendix  B. 

The  present  report  is  preliminary  in  that  the  results 
available  so  far  are  rather  fragmentary;  but  further  trial 
calculations  are  planned  and  will  be  described  in  subsequent 
reports. 

2.   The  Problem. 


Let  z,  r  be  cylindrical  coordinates.   It  is  assumed 
that  the  flow  is  symmetric  about  the  z  axis  and  that  the 
shape  of  the  shock  is  given  by  the  equation 

(1)     z  =  G{r)  =  Gq  +  Ggr^  +  G^r  +  ••'  • 

In  the  region  ahead  of  the  shock,  given  by  z  <  G(r) 
(to  the  left,  in  Figure  l) ,  the  flow  quantities  are 
constant  and  the  velocity  is  parallel  to  the  z-axis.   The 
air  is  assumed  to  behave  as  a  perfect  gas  with  the  equation 
of  state 
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(2)  P  =  (t  -  l)Sp> 

where  p,  &,      and  p,   are  respectively  the  pressure,  the 
specific  internal  energy,  and  the  density,  and  7  is  a 
constant.   It  is  convenient  to  take  the  units  of  length, 
density,  and  pressure,  (these  are  dimensionally  independent) 
to  be  the  radius  of  curvature  of  the  shock  at  the  axis, 
the  density  ahead,  and  the  pressure  ahead.   Then  the 
conditions  ahead  can  be  written  as : 

P  =  1, 

P  =  1, 

(5)  u  =  U, 

V  =  0. 

Here,  u  and  v  are  the  axial  and  radial  components  of 
fluid  velocity  and  U  is  a  constant  related  to  the  Mach 
number  M  by  the  equation: 

u  =  M/y-  . 

Behind  the  shock,  the  flow  variables  depend  on  z 
and  r.   It  is  convenient  to  introduce  coordinates  x  and 
y,   defined  by 


X  =  z  -  G(r),  z  =  X  +  g(y), 

(r) 

y=r-rQ,  r=y+r. 
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where  r   Is  a  constant  and  g(y)  =  G(y  +  r  ),   and  to 
expand  in  powers  of  x  and  y.   y  Is  a  radial  coordinate 
measured  from  the  point  of  expansion  (not  from  the  axis) 
and  X  Is  a  coordinate  measured  parallel  to  the  axis  from 
the  shock  (not  from  a  fixed  plane) .   The  shock  Is  given  by 
X  =  0.   If  f  stands  for  any  one  of  the  flow  variables 
u,  V,  p,  p,   we  write 

(5)     f  =  f(x,  y)  =y-       f  .  ,  x-^  y^ 

^  (J,  k)  ^   ^ 

Summed  over  all  pairs  of  non-negative  Integers   (j,  k)^ 
the  f  .  ,   are  the  coefficients  of  the  expansion.   This  Is 
understood  to  apply  in  the  region  x  >  0?   and  f (O,  y) 
Is  understood  to  denote  the  value  of  the  quantity  f 
Immediately  behind  the  shock. 

To  describe  the  Inclination  of  the  shock  at  point  y^ 
let  a  be  the  angle  between  the  shock  normal  and  the 
direction  of  the  Incident  flow.   Then 


sin  a  = 


s'(y). 


\/i  +  [g'(y)]^ 
(6) 


coe  a  = 


\/i  +  [g'(y)]^ 

Immediately  behind  the  shock,  the  flow  quantities  are 
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determined,  as  functions  of  y,   by  the  Ranklne-Hugonlot 
equations : 


(7) 


U  /  P  -  1 


/TTT^      \/  1  -  i/p 


(8)     p  =  ^P  "  ^   ,   where   6  =  ^-LA 

e  -  p  Y  -  1 


> 


(9)     {^  -u)   +   g'v  .  . /  (p  _  i)(i  _  i/p)   , 


(10)     (U  -  ^)S'  -  V   ^   0   , 
'l  +  g'^ 


where  g'   stands  for  g'(y)   and  u,  v,  p,  p,   stand  for 
u(0,  y),  v(0,  y),  p(0,  y) ,  p(0,  y) .   Equation   (7)   gives 
the  mass  flowing  through  unit  area  of  the  shock  in  unit 
time  (the  density  and  pressure  ahead  are  unity);   (8)   Is 
the  so-called  Hugonlot  curve.  In  the  pressure-density  graph, 
for  an  ideal  gas;   (9)   gives  the  Jump  of  the  normal 
component  of  fluid  velocity;  and   (lO)   states  the  conti- 
nuity of  the  tangential  component  of  fluid  velocity.   The 
irrationalities  are  easily  eliminated  from  equations   (7) 
to   (10):   from   (7)   by  squaring  both  sides,  from   (9)   by 
use  of   (7)>   and  from   (lO)   by  simply  cancelling  out  the 
denominator.   Therefore,  the  expansions  of  u,  v,  p,  p,   in 
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powers  of  y  can  be  obtained  by  these  equations  from  the 
expansion  of  g(y),   which  Is  assumed  known.   This  part  of 
the  calculation  follows  more  or  less  standard  lines  and 
will  not  be  described  In  detail.   As  result  of  It,  there 
Is  stored  in  the  machine  (actually  on  magnetic  tape) ,  a  dozen 
or  so  of  the  coefficients  of  each  of  the  flow  quantities: 
specifically,  the  values  of 

^0  k  '   ^0  k  '   PO  k  '   Pd  k  '     k  =  0,  1,  ...,  K. 

For  X  >  0,   the  quantities   u(x,  y) ,  v(x,  y), 
p(x,  y),   and  p(x,  y) ,   satisfy  the  partial  differential 

equations  t 

(11)  p[(u  -  vg')  ^  +  V  <^]  +  £^  =  0   , 

<)x     ^y  ()x 

(12)  p[{u   -   vs<)  ^  +  V  ^]   +  ^  -   g^   ^  =  0      , 

<)x     ^y    c)y      <)x 

(15)   -^  pu  +  (-^  -  g'  -^)pv  +   P^   =  0   , 
c)x       <)y      <)x      y  +  r^ 

(14)    'Yp[(u  -  va')°^  +  V  <^]  -  p[(u  -  vg')<^  +  V  ^]  =  0  . 

^x     <)y  ^x     ^y 

(ll)  and   (12)   are  the  equations  of  motion^   (Ij)   is  the 
equation  of  contlnuityi  and   (l4)   is  the  energy  equation, 
from  which  O  has  been  eliminated  by  use  of  the  equation 
of  state,   (2).   The  quantities  in  the  square  brackets,  [], 
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are  derivatives  along  the  streamlines;  the  terms  In  v  g' 
come  from  the  non-orthogonality  of  the  coordinate  system. 
As  before,   g'   stands  for  g'(y).   In  place  of   (1^), 
we  might  have  written  simply  that  the  derivative  of  pp"^ 
along  a  streajnllne  vanishes;  we  have  chosen  the  form   (1^) 
Instead  to  avoid  Irrational  operations. 

It  will  appear  that  there  Is  a  unique  power-series 
solution  of  the  equations   (11)   to   (l4)   that  satisfies 
the  conditions   (7)   to   (10)   for  x  =  0.   In  terms  of 
this  solution,  the  stream  function  tj/     is  given,  except  for 
a  factor  2ir,      by  the  equation: 

(15)   -^  =   T^(x,y)  =  I  U(y  +  r^f   -  (y  +  r^)  /  p(x}y)  v  (x}y)dx' 

The  streamline  t^  =  0  presumably  consists  of  two  parts:  the 
X-axis  and  a  curve  intersecting  the  x-axis  at  a  positive 
value  of  x;  the  second  part  may  be  taken  as  defining  the 
surface  of  the  body.   The  problem,  then,  is  to  find  this 
streamline  and  the  flow  between  it  and  the  shock. 

To  ascertain  the  region  of  the  variables  x,  y,   for 
which  the  series  converge  we  would  have  to  know  the 
singularities  of  u,  v,  p,  p,   regarded  as  analytic  functions 
of  two  complex  variables  x  and  y.   We  have  taken  the 
attitude  that  this  will  normally  be  beyond  one's  powers  of 
analysis  or  intuition  and  that  hence  one  will  have  to  be 
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satisfied  with  Indirect  empirical  evidence  of  convergence. 
The  nature  of  this  evidence  will  be  discussed  later,  but 
it  may  be  pointed  out  that  In  many  analogous  problems  In 
Incompressible  flow  the  region  of  convergence  Is  ample  to 
permit  calculating  the  flow  all  the  way  up  to  the  body  and 
even  to  extrapolate  It  analytically  for  a  considerable 
distance  Into  the  body.   An  example  Is  the  two-dimensional 
flow  obtained  by  superposing  a  uniform  flow  on  a  flow 
diverging  symmetrically  from  a  point  source.   In  this 
problem,  also,  the  streamline  i^  =  0  consists  of  two 
parts;  the  axis  (i.e.,  the  line  through  the  point  source  and 
parallel  to  the  flow)  and  a  curve  Intersecting  the  axis 
upstream  from  the  point  source.   The  second  part  may  be 
taken  as  the  surface  of  a  body  and  we  then  have  an 
Instance  of  an  incompressible  flow  past  a  two-dimensional 
blunt  body  or  rounded  wedge.   The  only  singularity  of  the 
complex  potential  is  at  the  point  source,  well  inside  the 
body.   This  circumstance  encouraged  us  to  believe  there 
may  be  an  adequate  region  of  convergence  in  the  detached- 
shock  problem  --  it  shows,  in  any  case,  that  the  stagnation 
point  (the  point  on  the  axis  at  which  u  =  v  =  O)   is  not 
generally  a  singularity. 

The  fear  has  been  expressed  that  there  may  be 
singularities  on  the  sonic  line,  which  separates  the 
subsonic  part  of  the  flow  from  the  supersonic  part.   Our 
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results  suggest  strongly  that  this  is  not  the  case,  at  least 
for  the  part  of  the  sonic  line  that  we  have  been  able  to 
explore  so  far.   Furthermore,  it  follows  from  equations 
(7)   to   (10)   that  the  expansions  In  y  for  x  =  0 
converge  well  beyond  the  sonic  point  on  the  shock,  in 
typical  cases. 

;$.   Method  of  calculation. 

The  procedure  to  be  described  will  determine  as 
many  as  one  wishes  of  the  coefficients  u.  ,  ,  v.  ,  ,  p.  ,  ,  p.  ,  , 
for  j  =  0,  1,  ...,  k  =  0,  1,  ...  .   Specifically,  given  a 
positive  integer  K,   the  coefficients  for  which  J  +  k  <  K 
will  be  determined.   In  one  calculation,   K  was  taken  to 
be  l4   (this  gave  120  coefficients  for  each  of  the  four 
flow  quantities);  later,   K  was  taken  as   10   (there  were 
then  66  coefficients  for  each  of  the  quantities). 

The  equations   (11 )   to   (l4)   are  first  solved 
for  the  x-derivatives ; 

,     P(u  -  vg.  )v[^  +  f^rr-]    H-  YP[pu^  +  g.^  -  pv^] 
^"^  YP(1  +  g'^)  -  p(u  -  vg')^ 

<)x     p(u  -  vg') 
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(^x        p(u  -  vg') 

(13)  JP   P[(u-vg-)i|.vif]-YPvi^ 
<)x  ypCu  -  vg') 


The  general  Idea  of  the  method  can  be  explained 
with  the  aid  of  Figure  2.   Each  dot  represents  a  pair  of 
values  of   J  and  k  for  which  we  wish  eventually  to 
know  the  coefficients  u.  ,  ,   etc.   Suppose  that  at  some 
stage  of  the  calculation  we  know  the  values  of  all 
coefficients  of  u,  v,  p,   and  p  whose  Indices  satisfy 
J  <  j'   for  some  j'   In  the  Interval   (0,  K)   as  well  as 
the  Inequality  J  +  k  <  K.   The  corresponding  dots  In 
Figure  2  lie  In  the  region  AAAA   (j*  =5,  for  the  case 
drawn).   Sor  brevity,  we  shall  say  that  the  functions 
u,  V,  p,  p,   are  known  In  that  region.   Such  a  region  has 
the  property  that  If  two  functions  are  known  In  It  their 
product  or  sum  or  difference  or  quotient  can  be  found  In 
the  same  region  by  formal  multiplication,  addition, 
subtraction,  or  division,  of  power  series.   Addition  and 
subtraction  are  trivial.   In  multiplication,  the  coefficient 
at  point  X  in  the  figure,  for  example,  of  the  product  of 
two  functions   f   and  g   involves  the  coefficients  of 
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f  and  of  g  at  points  lying  in  the  dotted  rectangle, 
of  which  the  point  X  and  the  origin  are  opposite  corners. 
A  similar  remark  holds  for  division,  which  is  performed  by 
the  rules  of  long  division  for  polynomial  expressions,  as 
taught  in  elementary  algebra.   Subroutines  have  been 
coded  for  performing  these  operations.   Further  subroutines 
are  available  for  formal  or  term- by- term  differentiation 
with  respect  to  y  and  integration  with  respect  to  x. 
These  subroutines,  and  the  manner  in  which  they  are  used, 
are  described  in  Appendix  B. 


A  B 


A  B 


Figure  2 
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If,  then,   u,  V,  p,   and  p,   are  known  In  the 
region  AAAA,   their  derivatives  with  respect  to  y  will 
be  known  in  the  shaded  region  of  the  drawing,  which  Is  the 
same  as  region  AAAA  but  shifted  downward  unit  distance 
because  differentiation  reduces  the  degree  of  each  term 
of  the  series o  By   use  of  the  operations  of  addition, 
multiplication,  and  division,  the  right-hand  side  of 
equation   (l6)   can  be  calculated  also  In  the  shaded 
region.   Similarly,  the  right-hand  sides  of  equations   (l?) > 
(18),   and   (19),   can  then  also  be  found.   At  this  stage, 
we  have  the  coefficients  of  the  four  partial  derivatives 
with  respect  to  x  in  the  shaded  regiono   Term-by-term 
integration  with  respect  to  x  then  gives  the  coefficients 
of  u,  V,  p,  p,   in  the  region   BBBB,   which  is  the  same 
as  the  shaded  region  but  shifted  unit  distance  to  the  righto 
If  we  now  merely  include  the  original  coefficients  of  the 
functions  on  the   k-axis,  we  have  information  of  the  same 
sort  as  presupposed  at  the  beginning  of  this  paragraph  but 
in  an  enlarged  region^  the  inequality  j  <  J'   being 
replaced  by  j  <  j'  +  1.   This  completes  what  we  shall  call 
a  cycle  of  the  calculation- 

The  information  available  at  the  beginning,  obtained 
from  the  Ranking- Hugonlot  equations .„  is  a  special  case  in 
which  the  region  kAJ\A     reduces  to  the  points  on  the  k-axis, 
corresponding  to   j '  =  0=   Therefore^  application  of  K 
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cycles  of  the  kind  described  above  completes  the  main  part 
of  the  calculation.   With  K  =  10,   about  2  hours  are 
required  on  the  UNIVAC 

After  the  main  part  of  the  calculation  has  been 
completed,  an  auxiliary  routine  Is  available  to  sum  the 
series,  calculate  certain  additional  quantities,  and 
print  the  results  In  edited  form  on  the  supervisory-control 
typewriter.   This  routine  automatically  repeats  Itself 
indefinitely.   At  the  beginning  of  each  repetition,  the 
machine  waits  for  the  operator  to  type  in  a  value  of  x 
and  a  value  of  y;  it  calculates  with  these  values;  it 
prints  the  results;  and  it  then  goes  back  to  the  beginning 
and  waits  for  another  pair  of  values  of  x  and  y  to  be 
supplied.   In  this  way,  the  operator  can  map  out  a  region 
of  the  flow  for  study.   The  quantities  printed  for  each 
point   (x,  y)   are  u,  v,  p,  p,  ■^,   and  three  check 
quantities.   The  stream  function  if     is  obtained  from 
equation   (15) I  the  integration  is  performed,  term  by  term, 
on  the  series,  before  summation. 

One  of  the  check  quantities,  called  the  Bernoulli, 
is  the  quantity 

(20)  B  =  -^:^  £  +  u^  +  v2  . 

Y  -  1  p 

It  follows  from  the  partial  differential  equations   (11 ) 
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to   (l4)   that  B  Is  constant  along  a  streamline.   (This 
Is  Bernoulli's  law.)   It  follows  from  the  Rankine-Hugonlot 
equations  that  B  Is  continuous  across  a  shock.   Clearly, 
B  has  the  constant  value 

(21)  B^^-^^1^ 

Y  -  1 

in  the  region  ahead  of  the  shock.   Therefore,   B  must 
equal  B   at  any  point  that  can  be  reached  by  a  streamline 
originating  ahead  of  the  shock.   This  applies  to  the  entire 
flow  up  to  the  surface  of  the  body  and,  by  analytic 
continuation,  to  the  extrapolation  of  the  functions  into 
the  body.   The  value  of  B  provides  a  sensitive  test  of 
the  errors  caused  by  truncation  of  the  series. 

The  other  two  check  quantities  have  to  do  with  the 
constancy  of  entropy  along  a  streamline,  where  p  should 
be  proportional  to  p'.   Prom  the  value  of  the  stream 
function  at  the  point   (x,  y),   the  machine  determines  the 
point   (0,  y-,)   at  which  the  streamline  through   (x,  y) 
crosses  the  shock.   Letting  p   and   p   denote  the  values 
of  pressure  and  density  at   (x,  y)  ,   and  p.,   and   p,   their 
values  at   (O,  y-,),   the  machine  calculates  all  these 
quantities  by  summing  the  series;  it  then  calculates  and 
prints  the  quantities 

(22)  p/p-^   and   (p/pj^)"^   , 
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which  should  be  equal.   This  provides  a  further  test  of 
the  adequacy  of  the  calculation. 

As  indicated  in  section  1,  the  series-summation 
subroutine  is  so  designed  that  it  omits  terms  whose  sig- 
nificance index  (niunber  of  significant  decimal  digits  in 
the  coefficient)  is  equal  to  0  or  1. 

4.   Results  of  the  trial  calculations. 

Near  the  very  beginning  of  the  calculation,  the 
machine  transfers  control  to  a  place  in  the  memory  where  it 
expects  to  find  a  routine  for  calculating  the  expansion  of 
g(y)   in  powers  of  y.   ^y  placing  a  suitable  routine  in 
this  portion  of  the  memory,  one  can  choose  the  shape  of 
the  shock  as  desired.   In  the  trial  problem,  the  shock 
surface  was  assumed  to  be  a  hyperboloid  of  revolution,  the 
asymptotes  having  just  the  right  slope  to  give  a  weak  shock 
(just  sonic)  at  large  distances  (as  expected  for  the  shock 
produced  by  a  finite  body) .   The  parameters  of  the  hyperbola 
were  so  chosen  that  the  radius  of  curvature  (of  the  shock, 
not  of  the  body)  was  unity  as  the  nose:  there  is  no  lack 
of  generality  in  this  choice  because  the  problem  scales. 

The  problem  was  further  characterized  by  the  choice 
of  Mach  number  M  =  12  and  the  ratio  of  specific  heats 
7  =  1.4  .   Several  calculations  were  made  for  this  problem, 
with  different  values  of  r  ,   to  obtain  series  expansions 
about  several  different  points  on  the  shock  surface. 
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The  location  of  the  streamline  -^  =  0  and  certain 
other  features  of  the  problem  are  Indicated  graphically  In 
Figure  5  (this  is  not  Intended  to  be  a  picture  of  the 
flowi  it  shows  more  curvature  than  a  picture  would,  because 
different  scales  have  been  used  on  the  two  coordinate  axes). 
Three  of  the  points  on  the  shock  about  which  expansions  were 
made  are  indicated  by  the  large  dots  labelled  "calc.  001", 
"calc.  002",  and  "calc.  003".    The  center  of  expansion 
for  the  fourth  calculation  was  at  y  =  0.8,   above  the  top 
of  the  graph.   The  shaded  bands  are  intended  to  indicate 
roughly  the  boundaries  of  the  regions  in  which  various 
degrees  of  accuracy  were  obtained.   In  the  region  between 
the  shock  and  the  left-hand  sl:^aded  band  the  flow  quantities 
are  accurate  to  around  six  decimal  places  or  betteri   as 
one  moves  away  from  the  shock  toward  the  body,  the  accuracy 
decreases.   The  estimation  of  accuracy  will  be  described 
briefly  belowi  it  is  of  course  highly  subjective  and  rather 
rough . 

In  Table  1  a   few  illustrative  results  are  listed. 
They  are  for  the  points   A,  B,  C,  D   in  Figure  J>,    which 
lie  on  a  horizontal  line  half  way  between  the  centers  of 
expansion  for  calculations   001  and  002s  the  quantities 
listed  in  the  table  were  computed  with  both  power  series 
for  each  of  the  four  points.   In  the  floating-decimal 
notation  used,  the  number  In  parentheses,  after  the 
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Figure  3 
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fractional  part,  Is  the  exponent  of  the  power  of  10|   the 
single  digit  after  the  exponent  Is  the  significance  Index  -- 
the  number  of  digits  of  the  fractional  part  that  are 
alleged  to  be  free  of  error  caused  by  normalization  In  the 
floating-point  subroutine. 

For  each  of  the  four  points ,  one  can  get  an 
Impression  of  the  accuracy  by 

a)  comparing  the  results  of  calc.  001  and 
calc.  002, 

b)  comparing  "Bernoulli"  with  the  exact  value 
208.6  from  equation  (21), 

c)  comparing  p/p-j^  with  (p/p-^)"^ , 

d)  observing  the  significance  Index. 

(In  connection  with  d),   It  should  be  noted  that  the 

number  of  correct  digits  cannot  exceed  the  significance  Index 

but  can  be  less.)   In  this  way  the  conclusion  was  reached 

that  most  of  the  quantities  listed  In  Table  1  are  correct  to 

about  7  decimal  places  for  points  A  and  B,   to  about 

5  places  for  point  C  and  to  about  2  places  for  point  D. 

In  Figure  4  the  pressure  p  Is  shown  graphically  as 
a  function  of  distance  x  from  the  shock  along  several  lines 
parallel  to  the  axis  at  different  radial  distances.   The 
pressure  distribution  and  other  features  of  the  flow  seem 
to  be  In  agreement  with  present  knowledge  of  such  flows. 

The  limitation  on  the  accuracy  In  the  present 
calculation  comes  apparently  from  the  limitation  Imposed  by 
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Figure  4 
loss  of  significance  on  the  number  of  terms  that  can  be 
used  in  the  power  series.   In  Table  2  the  significance 
indices  are  given  for  the  expansion  coefficients  of  the 
four  primary  flow  quantities,  for  calculation  001. 
It  is  seen  that  although  120  coefficients  were 
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calculated  for  each  quantity,  only  about  half  of  them 
could  be  used  during  the  series  summation.   In  later 
calculations  only  66  coefficients  were  calculated  for 
each  quantity.   The  problem  Is  being  recoded  for  double- 
precision  arithmetic  In  the  hope  that  this  will  permit 
retaining  more  terms  in  the  series. 

G.  B.  Whitham  has  suggested  that  by  altering  the 
formulation  of  the  problem  in  certain  ways  the  loss  of 
significance  can  possibly  be  much  reduced.   Other  choices 
of  dependent  and  Independent  variables  can  be  made;  and 
use  can  be  made  of  known  integrals  of  the  equations  so  as 
to  reduce  the  partial  differential  equations  to  a  third  or 
second  order  system.   If  this  is  done  without  introducing 
irrational  operations,  the  methods  described  above  can 
still  be  used  (only  the  pseudo-coding  portion  of  the  code 
has  to  be  changed,  and  this  is  a  relatively  small  task). 
Some  of  these  possibilities  are  now  being  explored. 
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APPENDIX  A  --  FLOATING- DECIMAL  SUBROUTINES 
WITH  SIGNIFICANCE  MONITORING 

Each  numerical  quantity  Is  represented  by  three 
numbers : 

1)  a  fraction  F,   which  lies  either  In  the 
Interval   -  1  <  F  <  -  1/10  or  In  the 
Interval  l/lO  <  F  <  1, 

2)  an  exponent   e,   which  is  an  Integer  In  the 
range  -  50  <  e  <  49, 

3)  a  significance  index  s,   which  is  a  non- 
negative  Integer  in  the  interval   0  <  s  <  9- 

The  limits  on   e  and  s   are  purely  arbitrary.   The  quantity 
represented  by  these  three  numbers  is  understood  to  be 

(A.l)  q  =  (F  +  10"®)  10^  . 

(According  to  this  definition,  the  last  significant  digit  is 
permitted  to  be  in  error  by  +1;   but  this  is  also  purely 
arbitrary:  for  example,  if  each  primary  datum  for  a  problem 
is  assigned  an  index  s  =  5  but  is  in  fact  correct  to  6 
decimal  places,  then  a  quantity   (F,  e,  s)   in  the 
calculation  will  be  correct  to  roughly  s  +  1  places,  if 
the  accumulation  of  round-off  and  other  sources  of  error  are 
negligible.)   Within  the  floating-decimal  routines,   F,  e, 
and  s  are  stored  in  separate  Univac  words,  but  elsewhere 
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In  the  problem  they  are  combined  or  composed  Into  a  single 
word.   The  composed  form  makes  use  of  the  twelve  characters 
of  a  Univac  word  as  follows;  the  first  9  of  the  characters 
are  the  first  9  digits  of  the  fraction   (F  +  l)/2   (this 
is  always  positive  and  hence  needs  no  sign) j  the  next   2 
characters  are  the  two  digits  of  the  integer  50  +  ej 
the  last  character  of  the  word  is  the  single  digit   So   Thus, 
the  numbers   1,  2,    -.01,   regarded  as  highly  accurate,  appear 
as   550000000519,  600000000519.   and  450000000^99,   re- 
spectively, while  TT,   rounded  off  to  four  decimal  places j, 
appears  as  657080000515- 

When  quantities  are  operated  on  arithmetically  by 
the  floating-decimal  routines,  the  significance  index  s   of 
the  result  is  determined  by  the  following  rules ; 

1)   Before  addition  or  subtraction  of  two  quantities, 
the  quantity  with  the  smaller  exponent  is 
denormalized  by  shifting  its  fractional  part  to 
the  right  a  number  of  places  equal  to  the 
difference  in  the  exponents  and  setting  its 
exponent  equal  to  that  of  the  other  number  (so 
far,  this  is  the  usual  floating-point  procedure)? 
at  the  same  time  its  significance  index  s   is 
also  Increased  by  an  amount  equal  to  the  number 
of  places  shifted,  or  up  to  the  value  3=9* 
whichever  is  the  smaller. 
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2)  As  the  next  step  In  addition  or  subtraction,  the 
significance  Index  of  the  result  Is  set  equal 

to  the  smaller  of  the  two  Indices  of  the  quantities. 

3)  If  addition  or  subtraction  results  In  overflow. 
In  consequence  of  which  the  resulting  F  Is 
shifted  one  place  to  the  right  and  the  value  of 
e  Increased  by  1,   the  value  of  s   Is  also 
Increased  by  1. 

4)  If  addition  or  subtraction  produces  cancellation, 
so  that  the  resulting  F  has  leading  zeros  after 
the  decimal  point,  the  quantity  is  normalized 

by  shifting  F  to  the  left  to  remove  the  leading 
zeros,  by  reducing  e  by  an  amount  equal  to  the 
number  of  places  shifted,  and  by  reducing  s  by 
the  same  amount  or  down  to  the  value  s  =  0, 
whichever  is  the  smaller.   If  s  becomes  zero 
In  this  process,   e  is  reduced  only  by  the  same 
amount  as   s  and  F  Is  set  equal  to  zero.   (The 
last  is  probably  not  essential,  but  prevents  apparent, 
though  false,  recovery  of  significance  by  accumulation 
of  debris  when  several  quantities  with  s  =  0  are 
added  together.) 

5)  In  multiplication  or  division,  the  significance 
Index  of  the  result  is  made  equal  to  the  smaller 
of  the  significance  indices  of  the  operands. 
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(The  single  shift  that  may  be  required  for 
normalization,  to  the  left  after  multiplication, 
to  the  right  after  division,  is  not  accompanied 
by  any  change  of  s) 

6)  Whenever  e  becomes  equal  to  -  50,   F  and  s 
are  both  set  equal  to  zero. 

7)  Whenever  e   exceeds   +  49,   the  machine  stops 
and  prints  "exponent  overflow"  on  the  supervisory 
control  printer. 

8)  Attempting  to  divide  by  a  number  that  has   s  =  0 
causes  the  machine  to  stop  and  print  "Illegal 
division"  on  the  supervisory-control  printer. 

If   s  =  0,   equation   (A.l)   shows  that  the 
information  on  q   is  so  vague  that   q  =  0  is 
one  of  the  values  consistent  with  this  information? 
division  by  zero  must  of  course  be  prohibited. 

When  the  subroutines  are  entered,  the  operands  are 
understood  to  be  in  the  Univac  registers  rP  and  rLj   in 
composed  form.   Two  additional  fictitious  registers,   P.O 
and  F.l,   each  consisting  of  three  memory  cells,  and 
containing  quantities  in  uncomposed  form  are  also  used.   The 
routines  are  entered  by  the  usual  Univac  instruction-pair 
R,  U|   the  routines  all  have  a  common  exit  (with  one  exception) 
but  there  are  several  entrances,  corresponding  to  the  various 
operations  that  can  be  performed.   The  possible  operations  are? 
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ADD:  (rP)  +  (rL)  to   rF  and  to  F«0 

SUBTRACTS  (rP)  -  (rL)  to   rP  and  to  P,0 

MULTIPLY;  (rP)  •  (rL)  to  rP  and  to  P.O 

DIVIDE;  (rP)  -f-  (rL)  to  rP  and  to  P.O 

MULTIPLY  AND  ADD;  (rP)  •  (rL)  +  (P.O)   to  P«0 

COMPOSE;  (P.O)   to  rP 

DECOMPOSES  (rP)   to   Pd 

SUBT  DECOMPD  NUMBER     (rP)  -  (P.O)  to   rP  and  to  P.O 

MULTIPLY  BY  INTEGER     (rP)  •   (rL)  to   rP  and  to  P.O 

DIVIDE  BY  INTEGER       (rP)  -r   (rL)  to   rP  and  to  P.O 

In  the  last  two  cases,  register  rL  contains  not  a  composed 
floating-decimal  n\imber  but  an  Integer  of  not  more  than  two 
digits   stored  in  normal  Univac  fashion  as  Integer  times 
10"   J  the  routine  converts  this  Integer  to  composed  form 
with  significance  index  s  =  9  before  performing  the  operations 
called  for. 

The  rules  given  above  make  no  attempt  to  allow  for  the 
accumulation  of  errors.   This  could  have  been  done,  following 
a  suggestion  attributed  to  L.  H.  Thomas,  by  carrying  along 
the  probable  error  of  each  quantity.  Instead  of  the  significance 
index,  and  always  computing  the  probable  error  of  a  result  on 
the  assumption  that  errors  of  the  operands  are  normally 
distributed.   This  would  have  been  slightly  longer,  and 
would  have  required  slightly  more  storage,  because  a  single 
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digit  cannot  give  a  sufficiently  accurate  measure  of  the 
probable  error;  It  was  not  deemed  to  be  worthwhile  for  the 
present  problem,  'where  the  main  danger  seems  to  come  from 
loss  of  significance  caused  by  cancellation. 

The  subroutines,  as  coded  for  the  Unlvac,  occupy 
151  words  of  storage  and  are  only  slightly  slower  than  other 
floating-decimal  routines  for  the  Unlvac. 

Similar  subroutines  are  being  coded  for  the  IBM  704^ 
but  using  double  precision.   The  effective  speed  of  the  704 
will  be  considerably  reduced  thereby;  but  one  may  hope  that 
the  reduction  will  not  be  too  severe,  because  some  use  can 
be  made  of  the  built-in  floating-point  facilities  of  the 
machine.   The  704  code  will  use  an  alternative  procedure 
suggested  by  N.  C  Metropolis;  Instead  of  keeping  a 
significance  Index  s,   one  supplies  each  floating-point 
number  with  Just  enough  leading  zeros  so  that  all  the  remaining 
digits  are  significant.   (More  likely,  the  digits  after  the 
leading  zeros  Include  a  fixed  number  of  guard  digits.,  at  the 
extreme  right  of  the  Txumher ,    to  guard  against  contamination 
of  the  truly  significant  digits  by  accximulatlon  of  rounding 
errors.   However,  this  Is  really  just  a  matter  of  point  of 
view  —  the  machine  has  no  way  of  knowing  which  digits  are 
considered  to  be  guard  digits  rather  than  significant  oneSo) 
The  Idea  of  Metropolis's  scheme  Is  that  there  should  be  no 
normalization  after  addition  or  subtraction  and  that  after 
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multiplication  or  division  the  result  should  be  so  shifted 
(in  the  case  of  multiplication,  the  product  is  supposed  to 
be  initially  in  a  double-length  shifting  register)  as  to 
have  exactly  the  same  number  of  leading  zeros  as  the  less 
accurate  of  the  two  operands . 

It  is  understood  that  the  GEORGE  machine,  which  is 
being  built  at  the  Argonne  National  Laboratory,  will  have 
a  form  of  significance  monitoring,  based  on  retention  of 
a  significance  index,  built  into  the  circuits. 

To  show  the  effectiveness  of  significance  monitoring, 
several  fragments  are  given,  in  Table  ^,    of  a  table  of 
sin  X,  computed  from  the  power  series  in  x  by  means  of 
the  Univac  floating-decimal  routines  referred  to  above. 
It  is  seen,  by  comparison  with  known  values  of  sin  x,   that 
in  all  cases  the  niomber  of  digits  claimed  to  be  significant 
is  rather  close  to  the  number  of  actually  correct  digits 
in  the  computed  values.   For  large  x,   as  is  well  known, 
some  of  the  individual  terms  of  the  power  series  are  very 
large,  although  the  sum  is  at  most  of  order  unity,  owing  to 
a  high  degree  of  cancellation.   For  angles  in  the  neighborhood 
of  1500  degrees,  the  situation  has  become  so  bad  that  values 
like   (0.  +  l.)-10   are  obtained  (note  that  the  indicated 
range  is  wide  enough  to  include  the  true  values).   Without 
significance  monitoring,  answers  like  0.xxxxxxxxx«10 
would  be  obtained,  where  the  x's   represent  digits  that  are 
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not  necessarily  zero.   What  happens  In  such  cases  is  that 
at  some  stage  of  the  summation  of  the  series  for  sin  x 
the  partial  sum  attains  a  value  of  the  order  of  10  -^j 
since  only  9  digits  are  carried  In  the  calculation, 
the  error  of  this  nixmber  Is  already  +  10  ,   and  subsequent 
calculation  will  not  in  general  reduce  the  error.   Such 
results  are  worse  than  useless:  it  was  noted  in  the  main 
part  of  this  report  that,  in  the  hydrodynamlcal  problem 
under  study,  the  coefficients  of  the  power  series  for 
u(x,  y),   for  example,  are  themselves  of  a  sufficiently 
complicated  nature  as  to  be  subject  to  serious  loss  of 
significance.   Therefore,  as  mentioned,  the  series  for 
u(x,  y)   behaves  like  an  asymptotic  series,  when  the  work  is 
performed  on  a  machine  with  a  finite  number  of  decimal 
places.   To  obtain  good  results  from  summing  the  series, 
it  is  essential  to  drop  systematically  those  terms,  and 
only  those  terms,  with  significance  below  a  certain  level. 
In  the  case  of  a  simple  calculation,  like  that  of 
sin  X,   there  are  many  well-known  remedies  for  loss  of 
significance.   In  that  Instance,  for  example,  one  can  replace 
the  angle  x   (as  is  usually  done)  by  an  angle  in  the  first 
quadrant   (0  <  x  <  Tr/2)   whose  sine  is  the  same  as  that  of 
X,   with  possible  exception  of  the  sign,  which  can  be 
supplied  separately:  the  series  then  converges  rapidly.   For 
many  of  the  large,  complicated  problems  for  which  the  modern 
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large,  fast  machines  are  suited,  it  Is  hard  to  see  how  one 
can  find  such  remedies  or  have  any  faith  in  results  obtained, 
without  the  aid  of  some  form  of  significance  monitoring 
whenever  floating-point  arithmetic  Is  used. 
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APPENDIX  B  --  SUBROUTINES  AND  PSEUDOCODE  FOR 
MANIPULATING  POWER  SERIES. 

The  basic  operations  with  power  series  were  described 
briefly  In  section  3  of  the  main  part  of  this  reports  they 
are  addition,  subtraction,  multiplication,  division, 
differentiation  with  respect  to  y,   and  integration  with 
respect  to  x',      they  are  all  purely  formal  operations 
performed  on  the  coefficients  of  the  series. 

If  a  function  f (x,  y)   is  represented  by  a  power 
series 


f(x,  y)  <-->Z:(j,  k)  fjk.J/   , 


where  J   and  k  are  non-negative  Integers,  we  define  the 
(K,  J ' ) -segment  of  the  series  as  the  set  of  coefficients 
f  .  ,   for  which   j  +  k  <  K  and   j  <  J'   (see  Figure  2).   In 
the  Univac  code,   K  is  restricted  to  be  not  greater  than 
l4.   With  that  restriction,  the  number  of  coefficients  in  a 
segment  cannot  exceed   120.   A  set  of   120  consecutive 
memory  locations  for  storing  the  coefficients  of  a  segment 
is  referred  to  as  a  pseudo-location.   The  code  recognizes 
three  pseudo-locations  in  the  Internal  memory  of  the  Univac  - 
these  are  called  pseudo-registers.   In  addition,  magnetic 
tapes  mounted  on  servos   3,  4,  ,..,  9  may  be  used  for 
storage  of  segments,  two  consecutive  blocks  of  6o  words 
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each  being  used  for  each  pseudo-location. 

To  represent  the  sum  or  difference  of  two  functions, 
corresponding  segments  are  added  or  subtracted  term  by  term, 
using  the  floating-decimal  subroutines  described  In  Appendix  A. 

To  represent  the  product  of  two  functions,  as 

f (x,  y)  g(x,  y)  =  h(x,  y) , 

the  corresponding  segments  are  combined  In  the  appropriate 
way,  viz. : 

'^  J  k  =  .f^Q     k-  =  0  ^  J'  k-  S  j  _  J.  k  -  k'  • 

This  can  be  a  rather  lengthy  calculatloni  therefore,  to  save 
time,  the  machine  takes  advantage  of  the  fact  that  the 
coefficients  h  .  ,   for  which  J  <  j '  -  1  have  already  been 
calculated  In  the  preceding  cycle  of  the  calculation  (see 
section  3) .   For  this  purpose  It  Is  assumed  that  the  order 
of  the  operations  on  the  power  series  Is  the  same  In  one 
cycle  as  In  the  next.   As  soon  as  the  segment  of  a  product 
has  been  computed.  It  Is  written  In  two  consecutive  blocks 
on  the  tape  on  a  servo  designated  as  servo  A  which  may  be 
either  servo  1  or  servo  2.   Similarly,  Just  before  calcu- 
lating a  product,  the  machine  reads  the  next  two  blocks  from 
a  servo  designated  as  servo  B  and  assumes  that  these  blocks 
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contain  the   (K,  j'  -  l)-segment  of  the  product,  so  that 
only  the  coefficients  with  j  =  j"   have  to  be  computed. 
After  completion  of  a  cycle,  these  two  tapes  are  rewound  and 
servos  A  and  B  are  automatically  interchanged. 

The  division  of  one  power-series  segment  by  another 
is  somewhat  more  complicated;  but  it  follows,  in  principle, 
the  long-division  algorithm  for  polynomials,  as  taught  in 
elementary  algebra.   The  presence  of  two  independent 
variables  complicates  the  work,  but  not  essentially.   Use 
is  made  here,  also,  of  information  obtained  in  the  preceding 
cycle.   It  is  necessary,  here,  to  carry  over  four  blocks  of 
information  (two  segments)  from  one  cycle  to  the  next,  by 
means  of  the  tapes  on  servos  A  and  B:  these  contain  the 
quotient,  as  developed  so  far,  and  the  remainder,  as 
developed  so  far.   It  is  also  necessary  to  bear  in  mind  that 
divisor  and  dividend  have  both  acquired  new  terms  since  the 
preceding  cycle,  so  that  the  first  step  in  the  division 
process  is  to  correct  the  remainder  so  as  to  take  account 
of  the  new  terms  in  the  divisor,  using  the  old  terms  in  the 
quotient.   The  usual  steps  of  the  long-division  algorithm 
then  follow.   It  is  necessary  to  assume  that  the  constant 
term  of  the  divisor  is   not  zero,  for  otherwise  one  would 
not  obtain  all  the  terms  of  a   (K,  j')-segment  upon 
dividing  one  such  segment  by  another.   In  consequence, 
expressions  such  as   (sin  x)/x  must  be  reduced  to  determinate 
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forms  before  the  coding  of  the  problem  begins. 

The  functioning  of  servos  A  and  B,  which  carry- 
Information  from  one  cycle  to  the  next,  need  not  be  the 
concern  of  the  programmer  (who  should  now  be  called  a 
pseudo-programmer) ,  as  they  are  taken  care  of  automatically 
by  the  subroutines. 

To  facilitate  the  use  of  the  subroutines,  an 
Interpretive  routine  Is  Included,  so  that  the  programmer 
(pseudo-programmer)  can  direct  the  manipulation  of  the 
power-series  segments  by  merely  writing  a  pseudo-code. 
This  Is  a  three-address  code  whose  operands  are  understood 
to  be  segments  of  power  series,  rather  than  numbers.   A 
pseudo-Instruction  normally  uses  eight  of  the  twelve 
characters  of  the  Unlvac  word  and  has  the  form   (OP  aa  PP  77) , 
where   "OP"   stands  for  the  operation  part,  which  may  be 
"AA"  for  addition,   "SS"  for  subtraction,  or  the  llkej  and 
where   "aa"  and   "pp"  normally  denote  the  addresses  of  the 
operands  and  "'yy"     normally  denotes  an  address  for  storing 
a  result.   Addresses   "01",  "02",   and   "03"   refer  to  the 
three  pseudo-registers  in  the  internal  memory;  addresses 
"30",  "40",  ...,  "90"   refer  to  the  tape  units  (servos)  of 
the  Unlvac*   Individual  blocks  on  the  tapes  are  not 

♦  "51",  "41",  ...,  "91"  are  used  when  writing  the  first  block 
on  a  tape.   This  causes  the  machine  to  read  a  sentinel  block, 
first  backwards,  then  forwards,  before  writing,  to  provide 
more  accurate  positioning  of  the  tape. 
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separately  addressed:  Instead,  the  tapes  are  always  written 
in  the  forward  direction  and  always  read  In  the  backward 
direction;  when  a  given  tape  unit  Is  addressed,  the 
address  automatically  refers  to  the  next  pseudo-location 
(next  two  blocks)  on  the  tape.  In  the  appropriate  direction. 
It  Is  up  to  the  pseudo -programmer  to  remember  which 
quantities  were  stored  In  the  successive  pseudo-locations 
on  a  tape. 

The  Instructions  of  the  pseudo-code  are  described  In 
Table  4,  where  they  are  classified  as  binary  operations, 
unary  operations,  and  red  tape.   Note  that  although 
differentiation  Is  a  unary  operation.  Integration  is  a 
binary  operation.   This  arises  from  the  necessity  of  getting 
the  constants  of  Integration  (the  coefficients  of  the  powers 
of  y  only)  from  somewhere  --  they  might  have  been  obtained 
from  the  preceding  cycle  of  the  calculation,  as  were  some 
of  the  coefficients  In  the  cases  of  multiplication  and 
division  but  are  In  fact  obtained  from  an  addressable 
memory  location.   The  special  multiplication  instructions  -- 
"multiply  by  constant"  and  "multiply  by  function  of  y 
only"  --  are  Introduced  merely  to  save  machine  time  by 
elimination  of  certain  steps  needed  only  when  two  general 
segments  are  multiplied  together.   In  each  case,  the 
invariable  part  (the  constant  or  the  function  of  y  only) 
is  automatically  taken  from  the  next  block  (only  one  is 
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needed  for  this  purpose)  of  the  tape  on  servo   "-"   (this 
Is  the  designation  of  the  tenth  servo  In  Unlvac  Instal- 
lations).  The  tape  on  servo   "- "  must  be  prepared  In 
advance  and  must  contain  the  constants  and  functions  of 
y   In  the  order  needed. 

The  subroutines  described  In  this  appendix,  together 
with  the  necessary  storage,  the  floating-decimal  routines, 
and  the  three  pseudo-registers,  occupy  about  920  words 
of  the  Internal  memory  of  the  Unlvac,  leaving  about  80 
for  the  pseudocode  Itself.   It  may  be  mentioned  that  the 
solution  of  the  partial  differential  equations   (l6)   to 
(19)   requires  just   71   words  of  pseudo-code;  their  use 
has,  of  course,  to  be  preceded  by  various  preliminary 
routines  (separate  memory  loads) ,  which  calculate  the 
coefficients  of  the  initial  segments  and  various  other 
functions  of  y,   and  prepare  the  ground  generally  for  the 
main  calculation. 
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TABLE  1    VALUES  OF  THE  PLOW  QUANTITIES 
AT  SELECTED  POINTS  ON  THE  LINE 
r  =  0.3  -  COMPARISON  OF  RESULTS  OP 
CALCULATION  001  AND  CALCULATION  002, 


QUANTITY 

CALC.   001 

CALC.   002 

X  =  .00  :  (POINT   A 

OF  FIGURE   3) 

u 

+.5^25005^7 

(+1)8 

+.342500557 

(+1)7 

V 

+.323105995 

(+1)7 

+.525105991 

(+1)7 

p 

+.153969787 

(+3)8 

+.155969795 

(+5)8 

p 

+.578120869 

(+1)9 

+.578120871 

(+1)9 

V' 

+.638936633 

(+0)9 

+.638956633 

(+0)9 

Bernoulli 

+.208600011 

(+5)8 

+.208600021 

(+5)8 

P/pl 

+.100000001 

(+1)8 

+.100000000 

(+1)8 

(p/pl)*^ 

+.100000001 

(+1)9 

+.100000001 

(+1)9 

X  =  .02  :  (POINT  B 

OP  FIGURE   5) 

u 

+.3017833^5 

(+1)7 

+.501785515 

(+1)7 

V 

+.319109025 

(+1)7 

+.519109157 

(+1)6 

P 

+.1569^7221 

(+5)8 

+.156947227 

(+5)8 

P 

+.580335245 

(+1)8 

+.580555255 

(+1)7 

if 

+.527254785 

(+0)8 

+.527254761 

(+0)8 

Bernoulli 

+.208600012 

(+5)8 

+.208600029 

(+5)7 

p/pl 

+.100462666 

(+1)8 

+.100462665 

(+1)8 

(p/pl)"^ 

+.100462665 

(+1)8 

+.100462665 

(+1)7 
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TABLE  1,    Cont'd 


QUANTITY 

CALC.   001 

CALC.   002 

X  =  .04  :  (POINT   C 

OF  FIGURE   3) 

u 

+.263767373 

(+1)7 

+.263772079 

(+1)7 

V 

+  .31475''^121 

(+1)6 

+.314759845 

(+1)5 

P 

+.158520611 

(+3)7 

+.158519681 

(+3)7 

P 

+.578736295 

(+1)8 

+.578733473 

(+1)7 

if 

+.416992527 

(+0)7 

+.416992149 

(+0)7 

Bernoulli 

+.208600072 

(+3)7 

+.208600491 

(+3)7 

p/pl 

+.100002552 

(+1)7 

+.100001954 

(+1)7 

(p/pl)'^ 

+.100002483 

(+1)7 

+.100001917 

(+1)7 

X  =  .08  :  (POINT  D 

OF  FIGURE   3) 

u 

+.191077903 

(+1)^ 

+.193800897 

(+1)4 

V 

+.306520089 

(+1)5 

+.306647263 

(+1)^ 

P 

+.158938105 

(+3)6 

+.158327999 

(+3)4 

P 

+.568773229 

(+1)5 

+.566961963 

(+1)5 

if 

+.203204727 

(+0)5 

+.203196357 

(+0)5 

Bernoulli 

+.208654675 

(+3)5 

+.208638915 

(+3)4 

P/Pl 

+.974128018 

(+0)6 

+.970387562 

(+0)4 

(p/pl)^ 

+.974630377 

(+0)5 

+.970287873 

(+0)5 
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TABLE  2   VALUES  OP   s   (NO.  SIGNIFICANT  DIGITS)  IN 

COEFFICIENTS  OF  POWER  SERIES  FOR  CALCULATION 
001.   WHERE  NOT  LISTED,  s  =  0 


6 

6 

6 

6 

6   1 

u(x,    y) 

6   2 

6   2 

8    3 

1 

6    3 

1 

1 

7    5 

2 

2 

1 

' 

6   4 

2 

2 

2 

1 

6   6 

4 

3 

2 

2    2 

6   4 

4 

4 

3 

2    2      2 

7   6 

4 

4 

4 

4    3     3    3 

8   7 

7 

5 

4 

4    4     3    3    3 

6 

6   1 

6 

6   1 

6   2 

V 

(x,  y) 

6   2 

1 

6    3 

1 

5   3 

2 

1 

1 

8   5 

2 

2 

6   4 

2 

3 

1 

1 

6   5 

3 

3 

2 

2 

6   3 

4 

2 

1 

2 

1 

1 

6   i^ 

5 

4 

3 

2 

2 

2 

7   5 

4 

5 

3 

3 

3 

2    1 

1 

7  6 

6 

4 

4 

4 

4 

5    3  3 

0  12 


->  J 


1 
2 

3 
3 
4 
4 

5 
5 
6 

5 
6 

8  7 


1 
1 
1 
2 
3 
3 
3 
3 
5 
5 
8 


2 
1 
2 
2 

5 
4 

5 


2 
2 
2 
4 
5 
4 


1 

3 

3 
5 
4 


p(x,  y) 


3 

2  3 

4  4 

^  3 


3 


1 
1 


1 

1 


1 
1 


1  1 


7 

7 

7 

8    1 

7    2 

1 

p(x,  y) 

7    3 

1 

7    2 

1 

7    4 

2 

7    3 

3 

2 

1 

7    5 

4 

2 

2 

7    5 

3 

2 

2 

1 

8    5 

4 

3 

2 

2   3 

7    5 

5 

3 

3 

3  2   2 

7   6 

5 

4 

4 

4  4   3     3 

9   7 

6 

5 

4 

4  4   3     3 
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TABLE  3 

TABLE  OF  SIN  x,    COMPUTED  3Y  POWER  SERIES,  USING  FLOATING- 
DECIMAL  ARITHMETIC  WITH  SIGNIFICANCE  MONITORING. 


X,  degrees 


sin  X 


00 
10 
20 
50 
40 
50 
60 
70 
80 
90 


+.100000000  (-49)  1 

+.173648179  {  00)  9 

+.342020143  (  00}  9 

+.500000003  (  00)  9 

+.642787611  (  00)  9 

+.766044445  {  00)  9 

+.866025405  (  00)  8 

+.939692621  {  00)  8 

+.984807759  (  00}  8 

+•999999991  (  00)  8 


500 
510 
520 
530 
540 
550 
560 


+.642793915  (  00)  5 

+.500000427  (  00)  5 

+.342001079  (  00)  5 

+.173656933  {  00}  4 

+.000000001  (-04)  0 

-.173646303  {  00)  4 

-.342015723  (  00)  4 


see  footnote 


1010 
1020 
1030 
1040 
1050 


-.887561759  {  00)  1 

-.882447235  {  00}  1 

-.797995997  00)  1 

-.699856927  (  00)  1 

-.538061153  (  00)  1 


1500 

1510 


+.000000001 
+.000000001 


03)  0 
03)  0 


*  When  the  significance  index  s   vanishes,  the  fractional 
part  F   is  set  equal  to  zero  in  the  uncomposed  form  of 
the  number  in  the  machine.   The   1   in  the  9  th  place 
in  the  printed  number  in  these  cases  results  spuriously 
from  composing  and  subsequently  decomposing;  it  could 
have  been  eliminated  by  a  slight  refinement  of  the  coding. 
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TABLE  4   PSEUDO- INSTRUCTIONS 

I,  Binary  operations   (OP  a  p  y) .   a,  p,  y  may  be  addresses 
of  pseudo-registers   ("01",  "02",  "0^")   or  of  tapes 
("50",  "40",  .«.,  "90").   Y  may  also  be   "31".  •••> 
or   "91"  for  writing  1—  block. 

A.  Preparation 

1)  contents  of  a  are  sent  to  01 

2)  "     "   P  are  sent  to  02 

B.  Main  stage 

(AA  -  -  -):  contents  of  01  plus  contents  of   02 

are  sent  to  01 5  contents  of  05  are  undisturbed. 

(SS  -  -  -) :  similar 

(MM  -  -  -):  1)  contents  of   01   are  sent  to   05j 
2)  contents  of  02   times  contents  of  05  to 
01 J   original  contents  of  01  remain  in  05? 
contents  of  02  are  undisturbed. 

(OD  _  _  _);  contents  of  01  become  the  numerator, 

contents  of  02  become  the  denominator?  division 
takes  placei  remainder  is  sent  to  05,   quotient 
to  Oil   contents  of  02  are  undisturbed. 

(Ix  -  -  -)%    integral  of  contents  of  01  with  respect 
to  X  is  sent  to  01 j  constants  of  integration 
are  taken  from  02 1   contents  of  02  and  of  05 
are  undisturbed. 
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TABLE  4,  Cont'd 

C.  Concluding  stage.   Contents  of   01   (result)  to  y. 

II.  Unary  operations   (OP  a  3  7)   (up  to  J>     per  pseudo- 
Instruction)  .   The  3  address  fields  are  associated 
with  pseudo-registers   01,  02,  05,   respectively.   A 
unary  operation  is  performed  involving  pseudo-register 
01  and  pseudo-location  a  if  and  only  if  the  address 
"a"   is   4=  0-   Then  similarly  for  02  and  p,   and  for 
03  and  Y,   in  this  order.   Unless  the  addresses 
involve  cross-references,  the  operations  are  independent. 
We  describe  below  the  cases  in  which  a,  P,  y     are  all  1=  0. 

(BB  a  p  y) :  Contents  of  a  are  sent  to  01 

II  II        Q        II        II        II        QO 

II  II        y  tt        II        II        Q-z 

(HH  037):  Contents  of  01  are  sent  to  a 

II         II      ry~)  11      M      II      Q 

II         tl      Q-2      II      II      II      y 

(Ey  a  3  7) :  v^  (contents  of  a)  is  sent  to  01 
II   /    II      "    ft^  "    "    "    02 

II      /         II  II  \    11        II        It        Q^ 

III.  Unary  multiplications.   (OP  a  —  y) .   Used  for  multipli- 
cation by  quantities  (constants  or  functions  of  y  only) 
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that  do  not  change  from  cycle  to  cycle.   The  address 
In  the   p  address  field  of  the  instruction  would 
normally  be  zero  -  it  is  in  any  case  ignored  by  the 
pseudo-machine. 
A.  Preparation; 

1)  Contents  of  a  are  sent  to  01 

2)  "     "  servo   "-"  are  sent  to  02 
(MC  a  0  y) I    Multiply  by  a  constant. 

(MP  a  0  y) ;     "     "   "  function  of  y. 

From  this  point  on,  these  operations  are  Identical 

with  MM  except  that  unnecessary  steps  are  omitted 

to  save  time. 
IV.  Red  tape 

(RW  -  -  -);  Rewind  servos   "— "   "1",  "2"?   and  inter- 
change servos   A  and   B   (A  =  "1",  B  =  "2"   or 
A  =  "2",  B  =  "1").   Servos   A  and   B  need  not 
concern  the  pseudo-programmer  except  that  initially 
their  tapes  must  be  filled  with  adequately  many 
blocks  of  zeros.   Initially,  servo   "— "  must  have 
the  needed  constants  and  functions  of  y  put  on 
the  tape  in  the  right  order. 

(CT  -  -  y)  t    Increase   J'j   if  then   j'  <  K,   pseudo 

machine  takes  next  pseudo-instruction  from  (Univac) 
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location  y;      If   J '  =  K,   Unlvac  stops,  but  then 
If  Unlvac  start  bar  Is  depressed,  Unlvac  takes  next 
Instruction  (not  pseudo-instruction)  from  Unlvac 
location  following  that  where  the  pseudo-instruction 
(CT  -   -   y)    Itself  is  located. 

(UU  -   -   y)  ''    Take  next  pseudo-instruction  from  Unlvac 
location  y. 

(RP  -  -  -):  Write  memory  dump  on  servo  B  (17  bl.) 
Advance  servo  A  17  blocks. 
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